function bcdf = custom_betacdf(p, a, b)
% custom_betacdf  - CDF of the beta distribution
%
% FORMAT:       bcdf = custom_betacdf(p, a, b)
%
% Input fields:
%
%       p           values for which the beta CDF is computed
%       a, b        beta function parameters
%
% Output fields:
%
%       bcdf        beta CDF of for (p, a, b)

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 3 || ...
   ~isnumeric(p) || ...
   ~isnumeric(a) || ...
   ~isnumeric(b) || ...
    isempty(a) || ...
    isempty(b)
    error( ...
        'BVQXtools:BadArgument', ...
        'Missing or invalid argument given.' ...
    );
end
if isempty(p)
    bcdf = [];
    return;
end
osize = size(p);
if numel(p) == 1
    if numel(a) > 1
        osize = size(a);
    elseif numel(b) > 1
        osize = size(b);
    end
end
if ~isa(p, 'double')
    p = double(p(:));
else
    p = p(:);
end
if ~isa(a, 'double')
    a = double(a(:));
else
    a = a(:);
end
if ~isa(b, 'double')
    b = double(b(:));
else
    b = b(:);
end

% build output
bcdf = zeros(osize);
bcdf(~(a > 0) | ~(b > 0) | isnan(p)) = NaN;
bcdf(p >= 1 & a > 0 & b > 0) = 1;

% regular cases
k = find(p > 0 & p < 1 & a > 0 & b > 0);
if ~isempty(k)
    if numel(a) == 1 && ...
        numel(b) == 1
        bcdf(k) = betainc(p(k), a, b);
    elseif numel(a) == 1
        bcdf(k) = betainc(p(k), a, b(k));
    elseif numel(b) == 1
        bcdf(k) = betainc(p(k), a(k), b);
    else
        bcdf(k) = betainc(p(k), a(k), b(k));
    end
end

% reshape
bcdf = reshape(bcdf, osize);
